function targetdiag(Bs,RMSDs,RMSDz,varargin)

% _________________________________________________________________________
% Call:
% targetdiag(Bs,RMSDs,RMSDz)
%
% Input:
%   Needed:
%   - Bs : Bias = <m> - <r> where m is model and r the reference field
%           or Normalized Bias: B* = uB = Bias / sig(r) where sig(r) is the
%           standard deviation of the reference field.
%
%   - RMSDs : signed unbiased Root-Mean-Square Difference
%               RMSD'=uRMSD=((1/N) * sum([(mi-<m>)-(ri-<r>)]^2) )^0.5
%                           * sign(stdm - stdr)
%             or RMSD*'=nuRMSD=Normalized RMSD'
%               nuRMSD= ( 1+ std* ^2 - 2 std* R ) ^0.5
%               where std* = stdm / stdr is the normalized standard
%               deviation, stdm the modelled standard deviation and stdr
%               the reference standard deviation.
%               R is the linear correlation coefficient:
%               R = ( (1/N) * sum[(mi-<m>)-(ri-<r>)] ) / stdm * stdr
%
%   - RMSDz : total Root-Mean-Square Difference
%             RMSDz = ( (1/N) * sum((mi-ri)^2) ) ^0.5
%               OR
%             other quantities if 'non_RMSDz' == 'on'
%
%
%   List of options:
%

%   - 'MarkerDisplayed' : 'colorar' (default) / 'marker'
%       -> marker: Experiments are represented by individual symbols
%       -> colorbar: Experiments are represented by a color described in a
%       colorbar
%
%   ------ OPTIONS considered only if 'MarkerDisplayed' == 'colorbar'
%   - 'non_RMSDz' : 'on' / 'off' (default)
%       -> Values in RMSDz are not corresponding to RMSDz but other
%       quantities
%   - 'title_RMSDz' : Title of the colorbar. Works only if 'non_RMSDz'='on'
%   ------
%   ------ OPTIONS considered only if 'MarkerDisplayed' == 'marker'
%   - 'name_experiment' : Name of the experiment used in the legend
%                 ([numbre of experiments x length of the longer name]
%                  array of characters)
%                       -- Or cell with experiment names
%                           >> name_experiment={'BATS','EUMELI'};
%   ------
%
%   - 'ticks' : define ticks position
%       -> default: 'ticks',[0:0.2:(maximum RMSDs)]
%   - 'subsampleticklabel' : draw a tick label for 1 over
%                            "subsampleticklabel" ticks
%       -> default: 'subsampleticklabel',4
%       Note: it is not used when ticklabelpos is defined.
%   - 'ticklabelpos' : position of the tick labels.
%       -> empty by default
%   - 'circles' : define the position of circles to draw (except circle = 1
%                 plotted by default when normalized diagram)
%       -> default : 'circles',[.7 (maximum RMSDs)]
%   - 'OU' : Observational Uncertainty
%       -> default : 'OU',0.
%   - 'normalized' : 'on'(default)/'off' -> Normalized target diagram or
%                    not.
%   - 'out' : 'print'/'screen' (default) -> Figure output
%   - 'filename' : EPS file name when printed output. Default file name is
%   'output.eps'
%
%
%
%   Example:
%   % --
%   targetdiag(Bs,RMSDs,RMSDz,'normalized','off')
%
%   References:
%   Jolliff et al., Summary diagrams for coupled hydrodynamic-ecosystem
%   model skill assessment, Journal of Marine Systems, 76 (2009) 64-82
%
%   taylordiag.m from G. Maze
%   (http://code.google.com/p/guillaumemaze/wiki//matlab_routines_taylordiag)
%
% ----------------------------------------------------------------------
% Created: G. Charria (08/2008)
% -- Updated: 03/2010 --
% v0 : Original code
% v1.0 : - option 'MarkerDisplayed'
% ----------------------------------------------------------------------
% _________________________________________________________________________

% __ Check inputs fields
narg = nargin - 3;
if mod(narg,2) ~=0
    error('targetdiag.m : Wrong number of input arguments')
end


% __ Define limits of the figure
maxx=max([RMSDs+0.2,1.2]);
maxy=max([Bs+0.2,1.2]);
minx=-max([RMSDs+0.2,1.2]);
miny=-max([Bs+0.2,1.2]);

maxi=max([maxx maxy]);
mini=min([minx miny]);

% __ Default
ticks=0:0.2:maxi;
subsampleticklabel=4;
%ticklabelpos=0.2:0.4:maxi;
ticklabelpos=[];
circles=[.7 maxi];
OU=0.;%OU=.8;
normalized='on';
out='screen';
filename='output.eps';
non_RMSDz='off';
title_RMSDz=[];
markerdisplayed='colorbar';
name_experiment=[];
overx=0;

for iopt = 1 : 2 : narg
    optname  = varargin{iopt};
    optvalue = varargin{iopt+1};
    switch lower(optname)
        case 'ticks'
            ticks=optvalue;
        case 'subsampleticklabel'
            subsampleticklabel=optvalue;
        case 'ticklabelpos'
            ticklabelpos=optvalue;
        case 'circles'
            circles=optvalue;
        case 'ou'
            OU=optvalue;
        case 'normalized'
            normalized=optvalue;
        case 'out'
            out=optvalue;
        case 'filename'
            filename=optvalue;
        case 'non_rmsdz'
            non_RMSDz=optvalue;
        case 'title_rmsdz'
            title_RMSDz=optvalue;
        case 'markerdisplayed'
            markerdisplayed=optvalue;
        case 'name_experiment'
            name_experiment=optvalue;
    end
end
close all
figure
if strcmp(out,'print')
    set(gcf,'Visible','off')
end
set(gcf,'Color','w')
hold on

% __ Remove default axes
set(gca,'XLim',[mini maxi],...
    'YLim',[mini maxi],...
    'XColor','w',...
    'YColor','w',...
    'XTick',[],...
    'YTick',[],...
    'XTickLabel',[],...
    'YTickLabel',[])

% __ Plot centered axes:
plot([mini maxi],[0 0],'k')
plot([0 0],[mini maxi],'k')

% __ Add circles
% _ 1 - reference circle if normalized
if strcmp(normalized,'on')
    theta=0:0.01:2*pi;
    rho(1:length(theta))=1;
    [X,Y] = pol2cart(theta,rho);
    %X=X+centerx; % centerx=0. = origin
    %Y=Y%centery; % centery=0. = origin
    plot(X,Y,'k','LineWidth',1.5);
end

% _ 2 - secondary circles
for it=1:length(circles)
    theta=0:0.01:2*pi;
    rho(1:length(theta))=circles(it);
    [X,Y] = pol2cart(theta,rho);
    plot(X,Y,':k');
end

% _ 3 - Observational Uncertainty (OU) threshold
theta=0:0.01:2*pi;
rho(1:length(theta))=OU;
[X,Y] = pol2cart(theta,rho);
plot(X,Y,'--b');

% __ Add ticks
for it=1:length(ticks)
    % X-Ticks
    plot([ticks(it) ticks(it)],[0 min([maxx maxy])/30],'k')
    plot([-ticks(it) -ticks(it)],[0 min([maxx maxy])/30],'k')
    % Y-Ticks
    plot([0 min([maxx maxy])/30],[ticks(it) ticks(it)],'k')
    plot([0 min([maxx maxy])/30],[-ticks(it) -ticks(it)],'k')
end

% __ Add ticks label
if ~isempty(ticklabelpos)
    for it=1:length(ticklabelpos)
        % X-Ticks
        text(ticklabelpos(it),-2*min([maxx maxy])/30,num2str(ticklabelpos(it)),'HorizontalAlignment','center',...
            'FontSize',8)
        text(-ticklabelpos(it),-2*min([maxx maxy])/30,num2str(-ticklabelpos(it)),'HorizontalAlignment','center',...
            'FontSize',8)
        % Y-Ticks
        text(-5*min([maxx maxy])/15,ticklabelpos(it),num2str(ticklabelpos(it)),'HorizontalAlignment','left',...
            'FontSize',8)
        text(-6*min([maxx maxy])/15,-ticklabelpos(it),num2str(-ticklabelpos(it)),'HorizontalAlignment','left',...
            'FontSize',8)
    end
    
else
    for it=1:subsampleticklabel:length(ticks)
        % X-Ticks
        text(ticks(it),-2*min([maxx maxy])/30,num2str(ticks(it)),'HorizontalAlignment','center',...
            'FontSize',8)
        text(-ticks(it),-2*min([maxx maxy])/30,num2str(-ticks(it)),'HorizontalAlignment','center',...
            'FontSize',8)
        % Y-Ticks
        text(-5*min([maxx maxy])/30,ticks(it),num2str(ticks(it)),'HorizontalAlignment','left',...
            'FontSize',8)
        text(-6*min([maxx maxy])/30,-ticks(it),num2str(-ticks(it)),'HorizontalAlignment','left',...
            'FontSize',8)
    end
end

% ___ Definition of markers
switch markerdisplayed
    case 'marker'
        kind=['+';'o';'x';'s';'d';'^';'v';'p';'h';'*'];
        colorm=['b';'r';'g';'c';'m';'y';'k'];
        n=1;
        marker(1:(size(colorm,1)*size(kind,1)),1:2)=' ';
        for ic=1:size(colorm,1)
            for ik=1:size(kind,1)
                marker(n,:)=[kind(ik,:) colorm(ic,:)];
                n=n+1;
            end
        end
        if (length(RMSDs) > 70)
            disp('You must introduce new markers to plot more than 70 cases.')
            disp('The ''marker'' character array need to be extended inside the code.')
            return
        end
        % __ Data points
        hp(1:length(RMSDs))=0.;
        for i1=1:length(RMSDs)
            hp(i1)=plot(RMSDs(i1),Bs(i1),marker(i1,:),'MarkerSize',8,...
                'MarkerFaceColor',marker(i1,2),...
                'Linewidth',2.5);
            if (strcmp(marker(i1,1),'o')||...
                    strcmp(marker(i1,1),'s')||...
                    strcmp(marker(i1,1),'d')||...
                    strcmp(marker(i1,1),'^')||...
                    strcmp(marker(i1,1),'v'))
                set(hp(i1),'MarkerEdgeColor','k')
            end
        end
        % __ Legend
        if ~isempty(name_experiment) && strcmp(markerdisplayed,'marker')
            
            a=get(gca,'Position');
            sfl=a(1)+a(3)+0.03;
            ll=legend(hp,name_experiment,'location','EastOutside');
            pleg=get(ll,'Position');
            overx=sfl+pleg(3)-1;
            if overx <= 0
                overx=0;
            end
            set(ll,'Position',[sfl-overx*sfl pleg(2) pleg(3) pleg(4)])
        end
            
    case 'colorbar'
        % __ Data points
        h=scatter(RMSDs,Bs,30,RMSDz);
        set(h,'MarkerFaceColor',get(h,'MarkerEdgeColor'),...
            'Marker','d')
        switch non_RMSDz
            case 'on'
                hc=colorbar('Location','NorthOutside');
                cp=get(hc,'Position');
                set(hc,'Position',[cp(1)+2*cp(3)/3 cp(2) cp(3)/3 cp(4)/2],...
                    'FontSize',8)
                if ~isempty(title_RMSDz)
                    title(hc,title_RMSDz)
                else
                    title(hc,'Color Scale')
                end
            case 'off'
                caxis([min(RMSDz) max(RMSDz)])
                hc=colorbar('Location','NorthOutside');
                cp=get(hc,'Position');
                set(hc,'Position',[cp(1)+2*cp(3)/3 cp(2) cp(3)/3 cp(4)/2],...
                    'XTick',[min(RMSDz) max(RMSDz)],'XTickLabel',{'Min. RMSD','Max. RMSD'},...
                    'FontSize',8)
                title(hc,'Color Scale')
        end
end

% _ to keep the ratio
axis equal;

% __ Axes Labels
ylabel('Unbiased RMSD','Color','k','Rotation',0,'Position',[maxi+3*maxi/30 0+3*maxi/30],'VerticalAlignment','middle')
xlabel('Bias','Color','k','Position',[0 maxi+3*maxi/30],'VerticalAlignment','middle')

% ------------------------------------
% ---- Resize step if needed
% ------------------------------------
fpos=get(gca,'Position');
fpos(3)=fpos(3)-overx*fpos(3);
fpos(4)=fpos(4)-overx*fpos(4);
set(gca,'Position',fpos)

% ------------------------------------
% __ Print eps file
% ------------------------------------
if strcmp(out,'print')
    set(gcf,'PaperOrientation','portrait')
    set(gcf,'PaperUnits','centimeters')
    set(gcf,'PaperType','A4')
    % Modification to display the whole legend    
    if strcmp(markerdisplayed,'marker')
        set(gcf,'PaperPosition',[1. 1. 38.+(10*(max((length(RMSDs)/10),3) - 3)) 38.+(10*(max((length(RMSDs)/10),3) - 3))])
    else
        set(gcf,'PaperPosition',[1. 1. 38. 38.])
    end
    %set(gcf,'PaperPosition',[1. 1. 18. 18.])
    print('-f','-depsc2','-r300',filename)
end
